#ifndef POINT_HPP
#define POINT_HPP

#include <limits>
#include <cmath>
#include <vector>


using std::sqrt;
using std::vector;


static const double
NA=std::numeric_limits<double>::infinity(); //should be quiet_NaN() but doesn't work

inline double min(double a, double b)
{
  return (a<b)?a:b;
}

inline double max(double a, double b)
{
  return (a>b)?a:b;
}

inline double min(double a, double b, double c)
{
  return min(min(a,b),c);
}

inline double max(double a, double b, double c)
{
  return max(max(a,b),c);
}

struct point //cartesian point with coupled euler's angles
  {
    point()
    {
      x=y=z=NA; //or may in would be numeric_limits<double>.quiet_NaN();
      o=a=t=0;
      
    }

    point(double _x, double _y, double _z,
	  double _o=0, double _a=0, double _t=0) 
    {
      x=_x;
      y=_y;
      z=_z;
      o=_o;
      a=_a;
      t=_t;
      
    }
    
    bool known() const
    {
      return (x!=NA);
    }
    operator float* () const
    {
     
      f[0]=x;
      f[1]=y;
      f[2]=z;
      return f;
    }

    operator double* () const
    {
      
      d[0]=x;
      d[1]=y;
      d[2]=z;
      return d;
    }

    double l2() const
    {
      return x*x+y*y+z*z;
    }
    
    double l() const
    {
      return sqrt(l2());
    }

    void normalize()
    {
      if(l2()!=1)
	{
	  double le=l();
	  x/=le;
	  y/=le;
	  z/=le;
	}
    }

    point norm() const
     
    {
      point p=*this;
      p.normalize();
      return p;
    }

    point operator-(const point &b) const
    {
      return point(x-b.x,
		   y-b.y,
		   z-b.z);
    }
    point operator+(const point &b) const
    {
      return point(x+b.x,
		   y+b.y,
		   z+b.z);
    }

    point operator*(const double &b) const
    {
      return point(x*b,
		   y*b,
		   z*b);
    }
    
    double x,y,z;

    double o,a,t;

    mutable float f[3];
    mutable double d[3];
    
    
    
};



class triang  //yet another triangle definition
{
public:
  point normal() const
  {
    if(size()!=3)
      throw;

    double x1=v[0].x;
    double x2=v[1].x;
    double x3=v[2].x;

    double y1=v[0].y;
    double y2=v[1].y;
    double y3=v[2].y;

    double z1=v[0].z;
    double z2=v[1].z;
    double z3=v[2].z;
    
    
    double
      vx1=x2-x1,
      vy1=y2-y1,
      vz1=z2-z1,
       
      vx2=x3-x2,
      vy2=y3-y2,
      vz2=z3-z2;
  

      
    double
      a=vy1*vz2-vy2*vz1,
      b=vz1*vx2-vx1*vz2,
      c=vx1*vy2-vx2*vy1;
  
    
    return point(a,b,c).norm();
  }
  
  point& operator[](int i){
    return v[i];
  }

  vector<point> v; 
  
};

  
typedef vector<triang> triar;

#endif
